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A theory of time-delayed coherent quantum feedback is developed. More specifically, we consider 
a quantum system coupled to a bosonic reservoir creating a unidirectional feedback loop. It is shown 
that the dynamics can be mapped onto a fictitious series of cascaded quantum systems, where the 
system is driven by past versions of itself. The derivation of this model relies on a tensor network 
representation of the system-reservoir time-propagator. For concreteness, this general theory is 
applied to a driven two-level atom scattering into a coherent feedback loop. We demonstrate how 
delay effects can qualitatively change the dynamics of the atom, and how quantum control can be 
implemented in the presence of time-delays. 
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Introduction. —Delayed autonomous feedback, where a 
signal is directly fed back to a system after a control¬ 
lable time-delay, is an important control tool for classi¬ 
cal systems m- It is highly attractive as a tool for 
stabilizing non-equilibrium states of fast dynamical sys¬ 
tems, where avoiding any time-costly signal-processing is 
crucial. Such stabilization is of great experimental and 
technological relevance HHS]. In particular, delayed au¬ 
tonomous feedback has been used to stabilize the high 
frequency dynamics of optical systems and high speed 
electrical circuits [Hill- 

Autonomous feedback is also receiving substantial and 
growing interest for controlling quantum systems HUS]. 
Because of the relatively short coherence time and fast 
dynamics of quantum systems, very fast feedback control 
possible with autonomous feedback is highly desirable. 
In addition, any measurement of the feedback signal will 
necessarily destroy its quantum character, making a fully 
quantum mechanical feedback loop that preserves coher¬ 
ence attractive from a fundamental point of view. Com¬ 
pelling evidence that this type of coherent feedback can 
outperform any measurement-based counterpart for im¬ 
portant quantum information processing tasks has been 
given HZIIIHI- 

A natural way of implementing coherent feedback con¬ 
trol loops is by coupling remote quantum systems via 
waveguides [T9H22] . Time-delays are unavoidable in prac¬ 
tice in such setups and are likely to become important if 
current experiments are scaled up to larger and more 
complex networks p51I25j . Despite of this, relatively lit¬ 
tle theoretical research has been done on delay effects for 
coherent quantum feedback. A major obstacle is the lack 
of tractable and general theoretical models for treating 
the highly non-Markovian dynamics induced by this type 
of feedback. The theoretical difficulty lies in the quan¬ 
tum correlations between the control target system and 
the in-loop quantum field: The field cannot simply be 
traced out, and one has to deal with a highly entangled 
quantum state over a continuum of degrees of freedom. 

Previous investigations have typically been limited to 


negligible delays [261127] , linear systems [IHl [28] , or sys¬ 
tems with special symmetries such as conservation of ex¬ 
citation number HSHSOj- For linear systems, some very 
promising theoretical demonstrations of the usefulness of 
delayed autonomous feedback to stabilize quantum sys¬ 
tems have been given recently. In [29] it was shown how 
it can be used to stabilize Rabi oscillations of an atom- 
cavity system in the single-excitation limit, and in |31j 
how it can enhance entanglement generated in a biexci- 
ton cascade in a quantum dot. Another study demon¬ 
strated that delayed coherent feedback might be used as 
a way of controlling the rate of convergence towards a 
non-equilibrium steady state in many-atom cavity quan¬ 
tum electrodynamics [5^ . 

In this letter we go beyond linear systems, and de¬ 
velop a general and tractable theoretical model for time- 
delayed coherent quantum feedback. This opens up re¬ 
search in a largely unexplored regime of quantum feed¬ 
back control. In particular, it allows for treating the 
important case of driven, non-linear systems, something 
which should be of immediate experimental relevance. 
We consider a generic setup where an arbitrary quantum 
system is coupled to a bosonic field forming a feedback 
loop. We show that the system’s density matrix can be 
found by evolving a time-propagator in an extended sys¬ 
tem space, followed by a generalized partial trace oper¬ 
ation. The evolution in this larger space is given by a 
differential equation for a time-propagator in Lindblad 
form. Interestingly, we can interpret this evolution as 
an unconventional quantum cascade |331134) , where the 
system is driven by past versions of itself. 

The derivation of our model uses so-called tensor net¬ 
work representations of quantum mechanical states and 
operators |5S|- These tools have their origin at the in¬ 
tersection of condensed matter and quantum informa¬ 
tion, where they are used to efficiently handle entangled 
many-body quantum systems. Recently, an intimate con¬ 
nection was made between continuum limits of certain 
tensor networks and output fields of open quantum sys¬ 
tems We develop these ideas further and find a 


2 


novel application of tensor networks in handling the dy¬ 
namics of a highly non-Markovian open quantum system. 
These developments could be of interest in themselves as 
a new approach to non-Markovian open systems theory. 

Below, we introduce the model putting emphasis on 
developing an intuitive picture of the dynamics. Techni¬ 
cal details are left to the Supplemental Material [35] • As 
a concrete example we consider a two-level atom coupled 
to a coherent feedback loop. We demonstrate two simple 
yet remarkable possibilities for delayed feedback control 
for this example: 1 ) spontaneous decay acting only for a 
controllable time, r, and 2) stabilizing Rabi oscillations 
far beyond the atoms coherence time in the absence of 
feedback. We discuss how these effects can be observed 
in a circuit quantum electrodynamics architecture m- 

Physical setup .—We consider a quantum system cou¬ 
pled to a single unidirectional bosonic field at two dif¬ 
ferent spatial positions, x = 0 and x = I, as depicted 
in |Fig. 1[ The field mediates a feedback loop for the 
system [IT]. We further assume that an arbitrary phase 
shift, (j), can be applied to the field between these two 
positions, such that the time-delay and phase are in¬ 
dependent parameters. The system-field Hamiltonian is 
H = Hs + Hb + V, where Hs is the system Hamiltonian, 
Hb = Jq dLJLJb^uJ)b(cJ) the free field Hamiltonian, and 
V the interaction Hamiltonian, 


V =i J duj^ ^ (ai 6 ^(w) — H.c.) 

+ i J (a 26 ^(a;)e-“"+*^ - H.c.), 


( 1 ) 


where r = l/c is the time-delay (c the speed of light), 
y/Ki ^2 is the coupling strength at the two positions, 
X = 0,1, respectively, oi and 02 are two system oper¬ 
ators, and H.c. stands for Hermitian conjugate. The 
field modes, 6 ( 0 ;), satisfy [ 6 (u;), 6 '^(w')] = S{uj — uj'). For 
generality, we allow the two system operators, Oi and 02 , 
to be different, but they could very well refer to the same 
operator—for example the dipole operator of a two-level 
atom or a cavity mode annihilation operator. The as¬ 


sumptions behind Eq. (1) are standard for open quantum 
systems, typically valid when the system is described by 
some frequency us ^ ki , 2 , see, e.g.. Ref. |i3| . 

To make the discussion more concrete, let us pause to 
consider a relevant example. A possible implementation 
is an optical cavity consisting of two mirrors, where the 
reflected field of one mirror is guided to be used as an 
input field on the other mirror (the inputs and outputs 
could be separated by circulators). In this case, one has 
the interaction in Eq. (1) with oi = 02 = a, for a system 
annihilation operator a, satisfying [a, a^] = 1 . ki _2 are 
in this example the linewidths of the two respective mir¬ 
rors. Hs describes the internal dynamics of the cavity, 
which could be non-linear due to the presence of other 
quantum degrees of freedom interacting with the cavity 



FIG. 1. (Color online.) Schematic of the setup. A nnidirec- 
tional bosonic field, E{x), interacts with the system, S, at 
positions a: = 0 and x = 1. The interaction at a; = 0 (* = 1) is 
with a system operator ai (aa) and a rate ki (fta). We assnme 
that an arbitrary phase shift, (f>, can be applied to the field 
between a: = 0 and x = I, snch that the time-delay, t = Ijc, 
and the phase are independently controllable. 


field. The equation of motion for the annihilation oper¬ 
ator in the Heisenberg picture can be found to be (see, 
e.g.. Ref. |42|): 


a(t) =i[Hs, a(t)] - - (ki -b K 2 ) a(t) - y/n^binit) 

I ( 2 ) 

- [^in(^ - t) + \f^\a{t - t)] . 

Here we have defined an input field &in(t) = 
-j= f°° due~'“^^*~*°^bo(w), where bo(u) are the initial 
values for b{u) in the Heisenberg picture. Eq. (2)| has 
the form of a delay differential equation [43], and makes 
the effect of the feedback quite clear. However, since the 
Heisenberg equations involve coupling between system 
and field operators, they are typically not efficiently solv¬ 
able in practice. Also, no corresponding master equation 
for the reduced system density matrix exists, in general, 
due to the finite time-delay. In the following, we present 
a practical scheme to integrate this type of dynamics by 
embedding the system in a larger space. 

A cascade of information from the past .—Our main 
result is a tractable model for the system dynamics af¬ 
ter eliminating the field degrees of freedom. The model 
suggests an intuitive picture where the system is driven 
by past versions of itself in a cascaded fashion. We here 
present the model and develop this picture, while leaving 
the technical details of the derivation to the Supplemen¬ 
tal Material [39] . 

To find the system state, ps{t), at an arbitrary time 
[k — I)t < t < kr, for k = 1,2,..., we evolve a time- 
propagator for a fictitious cascade of k identical copies 
of the system. The time-propagator, which we label 
Ss{t), is a superoperator on an extended system, 
and obeys a differential equation in the form of a cas¬ 
caded master equation as introduced by Carmichael and 
Gardiner [ 33 ] [ 33 ]. Note that we are here considering the 
master equation as a differential equation for the prop¬ 
agator, and not for a density matrix. Eor simplicity, we 
consider an incoming vacuum field, and we assume that 
the system and field are in a product state at time t = 0 
(the in-loop field is also initially in the vacuum state). As 
shown in the Supplemental Material [ 35 ] , the differential 
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FIG. 2. (Color online.) The time-propagator, £a{t), in Eq. (3) 
can be recognized as the propagator for a cascade of k iden¬ 
tical systems, Si. We can think of the copies as representing 
past versions of the system, i.e., the system is being driven 
by itself from the past. 


equation for the propagator then takes the form 

k 

(3) 

The integration variable, s, is an auxilliary time-variable, 
and the equation is to be integrated up to s = r, with 
the initial condition £o{t) = where I is the system 
identity super-operator. We have labeled k identical sys¬ 
tem copies by Si, I = l,...,fc. The superoperators H 
and V are defined by 


n[x],=[x ,.], 

V[X] •=x tX^ - ^x'fx • - ^ 
and the operators and T/,i+i are given by 


(4) 

(5) 


Him H.C.), (6) 


n.i+r-.^^s 

dO 


Li,i+i=y/Kia{’ + 


except for iJo,i — , j^k,k+i — 

y^e®‘^a 2 ^^ and Lk^k+i = The superscript de¬ 

notes the system on which an operator acts. Finally, 
we have defined A*^d(g) = for all I < k, and 

= 9[t — (k— l)r — where 9{s) is the Heav¬ 

iside step function, for any system operator A. 

1 is exactly the generator for a 


= H 


( 1 ) 


= H 


(fc) 


(7) 


Tn.i — 


The generator in Eq. 


cascaded chain of k identical quantum systems, as intro¬ 
duced by Carmichael and Gardiner [33] . An illustra¬ 
tion is given in Fig. 2 The evolution would describe a 
cascade in the usual sense if the time-propagator, £s(t), 
is applied to an initial state on the A:-fold system space, 
^(gifc However, the feedback problem is different, and the 
solution, psit), is found by imposing a peculiar type of 
“boundary conditions” on the propagator, as we will now 
explain. 

First of all, the integration variable, 


s, m 


Eq. (3) 


IS to 


be understood as a fictitious time-variable, and the equa¬ 
tion is to be integrated up to s = r, as already stated. 
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Eg. (9) 


3. (Golor online.) A diagrammatic illustration of 
Left: The propagator, £s{t), for the case fc = 3. The 


map is represented by a shape with lines attached to represent 
the input and output spaces. The labels indicate the systems 
associated to the lines. Middle: Eq. (9)| takes a state as input 
to system Si, while the output of system Si is mapped to 
the input of S 2 , and similarly the output of S 2 to the input 
of S 3 . The final output is a state for system S 3 . Right: For 
comparison, we show the application of the propagator to a 
state, psi ® PS 2 G) PS 3 , on the three-fold system space. This 
case corresponds to a conventional quantum cascade of three 
identical systems |33II34| (the choice of a product initial state 
is not essential). This diagrammatic notation is developed 
further in the Supplemental Material |39| . 


ps(t) is then found by acting with £r(t) on an initial state 
PSi (0) for the first system. Si , and essentially mapping 
the output of system Si to the input of system 5'i+i, for 
I = 1,... ,k — 1. The desired solution will be given as 
the output of system k, ps(t) = ps^(t). To explain this 
in more detail, we first have to introduce a generalized 
trace operation on the level of superoperators. For a su¬ 
peroperator, A, that acts on a tensor product of identical 
systems. S'! G • • • G 5^, we define the following generalized 
trace: 


Tr 


(Sc,Si) 


A* = 




(b|A( •(g)|b')(ii'l)b'z): (8) 


where \ii) and \iir) are bases for the two respective sys¬ 
tems, Si and Sii. Note that with I = V this is a partial 
trace, in the usual sense, but on the level of superopera¬ 
tors. More generally, this operation can be understood as 
mapping the output of system Si to the input of system 
Sp. 

We are now ready to write down an expression for 
ps{t)i given £r{t) found from Eq. (3)| 


psit) = Tr(5 ,_Sfc_i) • ■ ■ Tr(S 2 ,Si) ^T(t)pSi (0). (9) 


This equation, together with Eq. (3) constitute our main 


result, as they provide a practical scheme to find psit) for 
an arbitrary time t. In practice, the solution is thus found 
by first integrating Eq. (3) up to time s = t, and subse¬ 
quently computing the reduced state, psit), by acting on 
the initial state and taking the generalized partial trace 
. Eq. (9) To help build an understanding ofjEq. (9) 


we 


illustrate the trace operation diagramatically for the case 
fc = 3 in [Fig. 3| In the Supplemental Material, we give 
a derivation of Eqs. (§ and 0 using a tensor network 
representation of the time-propagator. 
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How can we now understand the dynamics induced by 
the feedback field? Eg. (3) [ suggest that the dynamics is 
given by a cascade of instances of the system, where each 
instance is driven by a past version of itself, from a time 
r earlier. What is highly non-trivial is that the feedback 
field that returns after a time r, is already quantum cor¬ 
related with the system it is driving. This leads us to 


Eg. (9) it is this equation that correctly account for the 
quantum correlations in time in the cascade picture. 

Delayed eoherent feedback for a two-level atom. —We 
illustrate the theory with a simple example: a two-level 
atom coupled to a coherent feedback loop. Both sponta¬ 
neous decay and resonance fluorescence through the feed¬ 
back loop is considered. This setup can, e.g., desribe an 
atom placed a (large) distance from a mirror, a problem 
with a long history in quantum optics (see Ref. [Hj and 
references therein). In the absence of a drive, the prob¬ 
lem can be solved analytically due to there being only a 
single conserved excitation between the system and the 
reservoir |45| . In the driven case, the problem has to the 
best of our knowledge previously only been considered in 
an approximate sense, employing perturbation theory in 
various limits [44) . 

The problem is defined by a system Hamiltonian Hs = 
£ (ct+ -I- (t_ ) , and coupling operators ai = 02 = U- . Here 
CT_ = \g){e\ is the atomic lowering operator, and a+ = 
(tT_)b £ is the drive amplitude, and we assume that the 
atom is driven on resonance for £ > Q. We take the rates 
to be identical, ki = K 2 = 7 , and assume a phase shift of 
(/) = TT in the feedback loop. 

Numerical results for the solution of Eqs. <§ and @ 
are shown in |Fig. 4[ The panels show three different 
cases: a) f /7 = 0, 6 ) f /7 = tt and c) f /7 = IOtt. The 
delay is chosen to be yr = 1.0 for the case £1^ = 0 , 
and otherwise equal to the Rabi oscillation period: t = 
27r/2£. The pink (light gray) lines show results with 
feedback, while the blue (dark gray) lines are analogous 
simulations without feedback, for comparison. 

We note two remarkable features in |Fig. 4[ First we 
consider the simplest case of spontaneous emission in 
panel a. In this case the atom decays exponentially to the 
ground state in the absence of feedback. In the presence 
of feedback, however, the feedback field starts driving the 
system after an initial transient period of time r, after 
which the population grows and eventually stabilizes at a 
steady state value. In steady state, destructive interfer¬ 
ence between two contributions to the output field, one 
coming from direct scattering and one from scattering via 
the feedback loop, prohibits the system from decaying. 
Hence, we have the possibility of letting the atom decay 
only for a controllable time. In steady state the system 
is dynamically decoupled from the decay channel. This 
phenomenon of feedback-induced dynamical decoupling 
of an atom from a decay channel has been demonstrated 
previously OlMlIlT]. 

Let us now look at non-zero drive strengths as shown in 



FIG. 4. (Color online.) Time-delayed coherent feedback con¬ 
trol of a two-level atom for three different parameter sets, 
as indicated in the figure. The pink (light gray) lines show 
the numerical results with feedback, while for comparison the 
blue (dark gray) lines show analogous simulations without 
feedback (ie., ki = 27 , K 2 = 0 ). 


panels b and c of Fig. 4| Here, the feedback induces long- 
lived Rabi oscillations, far beyond the coherence time of 
the atom in the absence of feedback. We have chosen r to 
coincide with the Rabi period, which is an optimal choice 
for stabilizing the Rabi oscillations. This means that t 
should be considered as a control parameter in its own 
right. In the bottom panel with S/y = IOtt and r = 0.1, 
the decay is extremely slow after the initial transient pe¬ 
riod of T. Numerical results have been verihed with a 
brute force numerical integration of the full system plus 
reservoir dynamics for small values of yr |48j . This was 
done by representing the feedback reservoir by a finite 
number of modes, truncated to have a small total pho¬ 
ton number. Such an approach however quickly becomes 
impractical for large yr (> 0.1). 

The simple example we have considered here could 
be realized experimentally in a variety of different plat¬ 
forms. A particularly appealing implementation is a 
circuit quantum electrodynamics architecture with an 
artihcial atom coupled to a one-dimensional waveg¬ 
uide [niiniiisiiso]. A meandering waveguide can be 
made to couple to the artificial atom at two different 
locations, or the artificial atom can be placed in a semi¬ 
infinite waveguide where the endpoint serves as a mirror. 
Such a setup was recently demonstrated experimentally 
in [^. A requirement to observe strong delay effects is 
yr > 0.1, which is readily achievable. In fact, signihcant 
delay effects are likely to be unavoidable even for mod¬ 
erate distances for strong coupling between the artihcial 
atom and waveguide. 

Conclusions .—We have shown that the problem of an 
arbitrary quantum system coupled to a coherent, held- 
mediated feedback loop, can be mapped onto a tractable 
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problem in a larger system space. This theory also yields 
an intuitive picture that helps to understand feedback 
mediated by a quantum field. For practical numerical 
integration, the approach presented here is superior to 
alternative approaches based on approximating the feed¬ 
back reservoir by a lower-dimensional system when the 
time-delay becomes comparable to the inverse linewidth 
of the emitting quantum system. 

The author thanks Howard Carmichael and Simon 
Whalen for many useful discussions and for performing a 
brute force numerical integration to verify the numerical 
results for small time-delays. The author thanks Alexan¬ 
dre Blais for helpful discussions. This work was sup¬ 
ported by NSERC. 
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This supplemental material to the paper “Time-delayed quantum feedback control” is organized 
as follows: First we give an introduction to a graphical “tensor network” notation for linear maps, 
based on Penrose’s tensor notation [1, 2], that is central to deriving the main results. We then 
briefly outline our approach to open quantum systems in a general context. Before approaching 
the feedback problem of the main letter, we consider the simpler problem of a system coupled to a 
conventional Markovian reservoir. We use this example to introduce the connection between open 
system dynamics and continuum limits of certain simple tensor networks. Finally, we consider the 
feedback problem of the main letter, and derive the main results stated there. 


I. A DIAGRAMMATIC DIRAC NOTATION 


In this section we introduce a graphical notation, which we can think of as a diagrammatic Dirac 
notation. Linear maps, such as vectors and operators, are denoted by shapes with “legs” attached 
to represent their input and output spaces. This type of notation was also developed at length in 
the context of open quantum systems in Ref. [2], and we refer to this work for more background 
information, although the focus there was on finite dimensional vector spaces. We here use slightly 
different conventions that are more convenient for our purposes. 

Let us first introduce some general notation. The Hilbert space of a quantum system A is denoted 
"Ha I where the subscript is used to differentiate between systems whenever there are more than 
one. Since we are dealing with open quantum systems, the state space of system A is a a subset of 
the Hilbert space of linear operators acting on PLa- We denote the latter L{'Ha), and refer to it as 
Liouville space. We can think of a state, p, either as an operator on PLa, or as a vector in Liouville 
space L{PLa)- To specify that a map, t, is an operator on PLa we write t : PLa —> PLa- To specify 
that it is a vector in Liouville space we write t 6 L{PLa)- It should otherwise be clear from context 
whether an object is thought of as an operator or as a vector. We therefore do not introduce 
any extra notation for the purpose of differentiating between an operator and the “vectorization” 
of that operator [12]. We also need to consider linear operators on L{PLa), which we refer to as 
superoperators. 

We first introduce the diagrammatic notation for an arbitrary complex separable Hilbert space, 
V, before specializing to the case V = L{PLa)- A vector u 6 V we denote by a shape with a leg 
going to the left: 



( 1 ) 


Similarly, a vector in the dual space w* E V*, we denote by a shape with the leg going to the right: 



( 2 ) 


Recall that the dual space, V*, is the set of all continuous linear functions, V —> C, and can be 
identified with V. An operator A on V we denote by a shape with two legs: 



(3) 


Finally, a number is denoted by a shape without any legs: 



( 4 ) 


^Electronic address: arne.loehre.grimsmo@usherbrooke.ca 
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We can equivalently think of the shapes as tensors, where each leg represents a tensor index, if we 
allow for countably infinite ranges for the indices [13]. 

The innerproduct on V, {w,v), is denoted 



(5) 


In general, composing objects is denoted by joining legs. Acting on a vector, v E V, with an 
operator. A, we denote 


Av = “EKD- (6) 

Acting on a dual vector would be denoted similarly, with the dual, say w*, placed to the left of 
A, and joining the two nearest legs. In tensor language, joining two legs corresponds to summing 
over a repeated index [2]. 

Tensor products of vectors, 6 W® V, are denoted by drawing the vectors vertically aligned, 
with one above the other: 


-© 

-© 


An arbitrary vector in the composite vector space, ri 6 W 0 V, is thus naturally drawn 


V - 
w- 


and an operator, C, on W 0 V: 


(7) 


( 8 ) 


C = 


V - 

w- 


>1 


- V 

w 


■W®V 


(9) 


In Eq. (8) and Eq. (9) labels are used to indicate the system assoicated to the legs. Such labes are 
typically omitted whenever they are clear from context. 

Finally, it is very useful to introduce a generalized trace operation for operators on tensor product 
spaces that allows us to connect an “output leg” to an “input leg.” Of course, connecting legs 
only makes sense if the systems associated to the legs are isomorphic. Consider the operator C in 
Eq. (9), and assume that V and W are isomorphic. We then introduce the following operation 


Tr(y.;,)C : A ^ T 
= '^x*Cyi, 


( 10 ) 


where X and y are any of the two spaces V or W, respectively. Xi and Ui are orthonormal bases 
for the two respective spaces, X and y. We note that this is a partial trace in the usual sense if 
X = y. We denote this operation diagramatically, by connecting the output leg for system X with 
the input leg for system y. For example: 


Tr 


(W.V) 


C: V^W 


( 11 ) 


In tensor language, this operation is of course nothing but summing over the corresponding repeated 
index. 

So far we have used a convention where objects are composed horizontally, right to left, and 
tensored vertically. We could equivalently have chosen to compose vertically, say top to bottom, 
and tensor horizontally. Below we take advantage of this freedom, and use different conventions 
for system and bath objects. 
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A. Density operators, superoperators and system-bath tensor products 


Let us now consider the case V = L{'Ha), for a quantum system A. Recall that the innerproduct 
on this space is given by {a a, ta) = We identify the operator-adjoint, with the dual 

vector a\ e L{'Ha)* of aA € L{'Ha)- 

It is convenient to introduce a special notation for the identity operator, I a, on the system 
Hilbert space 'Ha, viewed as a vector in Liouville space: 


Ia — —O G L{Ha), 

and its dual 

I*A = o- e L(Ha)*■ 

This allows us to draw for the trace of ta € L('Ha)’- 

trTA = {Ia,ta) = 


( 12 ) 

(13) 

(14) 


Later, when dealing with system-bath dynamics, it will be useful to have a notation that dis¬ 
tinguishes between system and bath objects. We now introduce a set of specialized conventions 
for this purpose. Consider a system-bath tensor product. S' ® B. System states are denoted as in 
Eq. (1): 

PS =HT-ts)- (15) 

We distinguish this from states of the bath, B, by drawing the latter with the leg going downwards: 

(^B L{'Hb), (16) 


A superoperator, Asb, that acts on the composite system-bath space we then draw 


B 

_ 1 _ 


Asb = s - Asb - s , 


(17) 


B 


where the labels indicate the system associated to the legs. This convention is practical for dealing 
with superoperators acting on product states, ps ® ^b, which we can now draw 


(18) 


AsbPs ® — Asb ~(p^- 


We need to compose such superoperators in two different ways. Superoperators, say Asb and 
Bsb', that act on two different bath systems, B and B': 


B' 

_X_ 


B 

I 


Bsb'Asb = s - Bs,b' — As,b - s . 


(19) 


B' 


\ 

B 


And, superoperators that act on the same bath system, B, but two different systems S and S’: 

B 


Bs'bAsb — 



( 20 ) 
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Finally, we consider how the generalized trace defined in Eq. (10) acts on the level of superop¬ 
erators. Consider a superoperator, A, acting on two isomorphic systems S and S': 


We then have that 


^ : L{ns) ® L{ns') ^ LiHs) ® L{ns') 


s ■ 

s'- 


.4 


■ s 

■s' 


Tr(s',5)^» : LiHs) ^ 




•®Ns')(is'l)|js), 


( 21 ) 


( 22 ) 


where {|*s)} s-nd {|is')}i orthonormal bases for the two respective systems. 


II. OPEN QUANTUM SYSTEMS 

In general, the evolution of an open quantum system is given by a Hamiltonian for the system 
and a bath, H{t) = Hs(t) + Hsit) + V{t), together with an initial state po for the composite 
system, usually assumed to be in product form po = Ps®Pb- Here Hs{t) is a system-Hamiltonian, 
HB{t) is the Hamiltonian for the bath, and V{t) is an interaction Hamiltonian. In general, we 
allow these Hamiltonians to be time-dependent. The reduced state of the system, ps{t), at a time 
t, can then formally be written 

ps(t) = Ub W(t)po, (23) 

where tra denotes a partial trace over the bath, and we have introduced a time-evolution super¬ 
operator defined through 

W(t)p = Texp i J pTexp J (24) 

Here T denotes the time-ordering operator. The method we introduce in the following sections is 
based on “trotterizing” W(t). That is, using a Suzuki-Trotter decomposition [3], we seek to write 
W(t) as a limit, 

W(t) = lim Un ■ ■-Ui, (25) 

N—^CiO 

where if„p = C/„pC/f is a unitary conjugation. For a common class of system-bath Hamiltonians, 
we show, in a sense that becomes clear below, that can be made arbitrarily close to the identity 
map as N is increased. Eq. (25) can thus be understood as a product integral [4] with a well-defined 
continuum limit as A —> oo. 

For finite N, the product in Eq. (25) can be represented diagrammatically using the notation 
introduced in the previous section. This can in certain cases expose a structure to the problem that 
would not easily be seen otherwise. This is indeed the case for the feedback problem considered in 
the main letter, and we exploit this to find a product formula for a time-propagator after tracing 
out the bath degrees of freedom: 

S(t) = lim (26) 

N—¥-oo 

where £„ is a superoperator on an extended system. S', that approaches the identity map as 
N ^ oo. The state ps{t) is finally found by acting with £{t) on an initial state ps(0) and 
performing a generalized partial trace, of the type introduced in Eq. (22). 

Although this approach is particularly useful for the non-Markovian feedback problem, it is 
instructive to apply these ideas first to the simpler problem of a Markovian bath, i.e., a conventional 
reservoir without any feedback loops. This is done in the next section, before we subsequently 
attack the feedback problem of the main letter. 
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III. MARKOVIAN DYNAMICS 

We assume that the system’s internal dynamics is described by a Hamiltonian, Hs, that we leave 
arbitrary. The coupling to the bath is assnmed to be linear {e.g., a dipole conpled to the electric 
field). More specifically, we assume a system-bath Hamiltonian, H = Hs + Hb + V, where (in 
units where ?i = 1) 

Hb = [ dcoujb^ {uj)b{uj) (27) 


is the bath Hamiltonian, and V is the interaction Hamiltonian, given by 


V 



— o^&(tj)) . 


(28) 


Here, k is a rate describing the coupling strength of the system to the field, a is an arbitrary system 
operator, and the field modes, b{uj), satisfy [&(w), &^(a;')] = 6{uj — uj'). 

It is now convenient to go to a rotating frame with respect to the bath Hamiltonian Hb- The 
interaction Hamiltonian then takes the from 


V{t)=i J (^ab^ — a^b{uj)e *“*) = z (L6^(f) — I/^6(t)) , (29) 

where in the last equality, we have introduced the system operator L = ^/na, and the time-domain 
field 


b{t) = ^l d^e—*6(cc), (30) 

which satisfies [b{t),b'^{t')] = 5{t — t'). Formally, the solution for the system state at a time t can 
be written 

ps(t) = trB>V(t)ps(0), (31) 

where the evolution superoperator VV’(t) is given through 

VV’(t)p = Texp i j dt'H'(t')^ pTexp J dt'H{t')^ , (32) 

where H{t) = Hs + V{t) is the Hamiltonian in the rotating frame with respect to Hb- 

The next step is to go to discretized time. This allows us to expose the structure of the evolution 
operator in terms of a “tensor network,” using the diagrammatic notation from Sect. I. We do this 
by first introducing ladder operators 


b 


n 


1 

\/Ai 




&(s)ds. 


(33) 


which satisfy the usual commutation relation [&„, b]^] = Snm- Here, = nAt, for n = 1,..., N, 
where t = NAt- Using the Suzuki-Trotter formula [3], the evolution operator can be written 


W{t) = lim UN---U 2 U 1 - 

N—yoo 


(34) 


Here, and in the following, V — > 00 refers to the continuum limit where NAt = t is held fixed. 
The unitary superoperator Un is defined through Un • = U„ • with 

Un = exp (^—iAtHs + \rAt{Lbl^ — L^b^j - (35) 

We now work with finite N, and eventually take the continuum, N —> 00 , limit. Eq. (34) and 
Eq. (35) then suggest the following picture: The system interacts with a collection of harmonic 
oscillators, labelled by n, where the n’th oscillator interacts with the system at time for a short 
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* * * > B >B >B >B >B • • • 


S 


FIG. 1: A “conveyor belt” of harmonic oscillators, labeled by n, moves past the system as shown in the 
figure. The n’th oscillator interacts with the system for a short time At. 


time At, before leaving the system. After this it never interacts with the system again. We refer 
to this as the “conveyor belt picture,” and it is illustrated pictorially in Fig. 1. 

The superoperator hl„ we denote diagramatically as 


71 

L{„= s 

-T 

71 



(36) 


where the labels indicate the system associated to the legs: S for the system and n for the n’th 
oscillator. The result of acting with this superoperator on a product state ps ® is then denoted 


UnPs ® cr„ = 


(37) 



We are now in a position to represent diagrammatically the system-bath state at time tjv, assuming 
an initial state ps ® cti ® (T 2 ® ® cfn' 


PSBitN) 




(38) 


As noted in the previous section, the shapes appearing in the diagram above, Eq. (38), can be 
thought of as tensors. A diagram consisting of many tensors, connected by joining their legs, like in 
Eq. (38), is therefore often called a “tensor network” [5]. In fact, the tensor network representation 
of the state in Eq. (38) is analogous to a matrix product state (MPS) [6], but with superoperators 
and density operators in place of operators and kets, respectively. MPS are famous for their role 
as a variational ansatz in the density matrix renormalization group method [6], used to find the 
low-energy spectrum of many-body Hamiltonians. In light of this similarity, we call a state on the 
form of Eq. (38) a superoperator product state (SPS) [14]. If taking the continuum limit (more on 
this below), N —> oo, the resulting state we can similarly call a continuous SPS (cSPS), in analogy 
with the recently introduces cMPS [7]. The connection between output fields of open quantum 
systems and cMPS has been exploited to use open quantum systems to generate cMPS states, in 
the context of quantum simulation, both in theory [8, 9] and in experiment [10]. 

In the present context, however, we are chiefly interested in the reduced system dynamics, 
Ps{t) = BPS sit)- That is, we want to trace over the bath degrees of freedom in Eq. (38). This 
gives a reduced system state (cf. Eq. (14)): 


Ps{tN) 




(39) 


We would next like to retrieve the continuum, Af —> oo, limit, to find a more compact form of 
Eq. (31) for ps{t)- To keep the discussion as simple as possible, we assume an incoming vacuum 
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field, i.e., we assume (7„ = |0)(0|, for all n. We can then find an analytical expression for the map 
f 

*-"n ? 



= exp {—iAtUlHs] • +AtV[L] •), 


where we have defined superoperators 

n[x].=[x ,.], 
v[x ]. =w • wt - 


(40) 

(41) 

(42) 

(43) 


In Eq. (40) higher order terms in At have been neglected. In the continuum limit, N —> oo, we 
define a time-propagator, £{t), which can be written as a product integral [4]: 


N 

£{t)=lim (44) 

N—^oo 

n=l 

For the case of a vacuum bath, this is simply 

£{t) =eyip {-im[Hs]+W[L]). (45) 

The Lindblad equation is retrieved by differentiating: 

(t) = -in[Hs]S{t) + V[L]£it). (46) 

The initial condition for the time-propagator is £^(0) = Is, where Is = Is * Is is the system 
identity superoperator. Finally, the system state, ps{t), is found through 


Ps(f) = £{t)ps{0)- 


(47) 


IV. A SOLUTION TO THE FEEDBACK PROBLEM 
A. A diagrammatic representation in discretized time 


Equipped with the graphical notation, and having warmed up with a Markovian reservoir in the 
previous section, we are ready to attack the feedback problem of the main letter. 

Recall that we are now considering the following interaction Hamiltonian: 


U y duj^ (ai6^(tj) — H.c.) 

+ i duj^ (a 2 &t(w)e-*“"+'^ - H.c.) , 


(48) 


where t = l/c is the time-delay, is the coupling strength at the two positions, x = 0,1, 

respectively, and ai and 02 are two system operators. We follow the steps of the previous section, 
and go to a rotating frame with respect to the bath Hamiltonian Hb = f d^JUJb^uJ)b(cAJ): 


V(t) =i J dtj^p-(ai5^(a;)e*‘^* — H.c.) 

+ i j dLo^J~^ - H.c.) 

= i (^Lib\t) - L\b{t)j + i (^L2b\t - t) - L\b{t - t)) , 


(49) 


( 50 ) 




f 



[ UaTH— ^ 


n — M 


FIG. 2: A “conveyor belt” of harmonic oscillators, labeled by n, moves past the system. The n’th and the 
(n — M)’th oscillator interacts with the system simultaneously, for a short time. At. 


where we have introduced the time-domain field, b{t), defined as in Eq. (30), and Li = y/Kiai, 
1/2 = e®'^yre 2 a 2 . The evolution superoperator, W{t) defined as before in Eq. (32), can again be 
written in terms of a product integral using the Suzuki-Trotter formula: 

W{t) = lim UnVn-m ■ ■ ■bl 2 V 2 -MUiVi-M- (51) 

N—^OO ^ ' 

Here = f/t and V„ = K • are unitary conjugations with 

Vn = exp -t i'/At(Libl 

Vn =exp -b i^/M{L 2 bl 

where is defined as in Eq. (33). NAt = t and MAt = t are both to be held fixed as A" —t oo. 

Eq. (51) and Eq. (52) suggest the following picture: The system interacts with a collection of 
harmonic oscillators, labelled by n, where the n’th oscillator interacts with the system at time 
for a short time At, before leaving the system. The same oscillator comes back to interact with 
the system a time t = AIAt later. After this it never interacts with the system again. This is 
illustrated pictorially in Fig. 2. 

We assume that any oscillator that has not yet interacted with the system is in a product 
state with everything else. We also assume that at time t = 0, the initial state is on the form 
Ps ® cTi-M ® ® fii (g • • • (g (jjv The evolution up to time r = Im = MAt is then essentially 

entirely analogous to what we had in the previous section. Diagrammatically, we can write 



PssitM) 



(53) 


where we have suppressed the labels on the oscillator states, an, for notational convenience. The 
oscillators with n < 1 will at this point never interact with the system again, and can therefore be 
traced out. Diagramatically this is represented by 


PSLitm) = 



(54) 


Here L stands for “loop,” and refers to the subset of bath-oscillators that have interacted with the 
system exactly once. Now, to integrate to time Im+i we need to attach the two superoperators Vi 
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andZ^M+i to Eq. (54): 


PSL{tM+l) 



(55) 


Notice that since Vi acts on the same oscillator as Ui, the upper (input) leg of the former has to 
be connected to the lower (output) leg of the latter. We can equivalently redraw this diagram in 
the following way: 


PSL{tM+l) 



(56) 


It is stressed that Eq. (55) and Eq. (56) are equivalent diagrams, drawn in two different ways. The 
latter, Eq. (56), is the preferred way of drawing the diagram, as will become clear in the following. 
In fact, this is the crux of the problem, and here the power of the diagrammatic notation comes 
into play. The diagrammatic notation is equivalent to the usual algebraic Dirac notation, but 
the same expression can be drawn in several ways, and the diagrammatic notation can therefore 
be richer than the algebraic notation. In some situations, as is the case here, the diagrammatic 
notation can reveal a structure to the problem that would not easily be seen otherwise. 

As we continue evolving the state for times t < tjsi < 2t, we keep attaching alternating V’s and 
W’s to the second row in Eq. (56). So, for example, the state at time t 2 M = 2t, we draw 


PSL{t2M) 



(57) 


At this point we can do the same trick as we did at time Im, and start a new row to attach the 
next pair of superoperators. That is, the state at time t 2 M+i we can draw 


PSL{t2M+l) 



(58) 


And so it continues: We evolve the state by attaching V’s and U’s, and start a new row at each 
n = kM + 1, where k = 1,2,.... Now, the system state psitn) = ti'LPsL(ijv) at time fw, is found 
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by tracing over the field. For example, the system state at time t^M+i is found from Eq. (58): 


Ps(t2M+l) 



(59) 


Notice that there is only one unpaired leg, going to the left, as it should be for a reduced system 
state. 

We are now ready to represent diagramatically a formal solution for the system state at an 
arbitrary time, (fc — l)r < tvr < kr, where fc = 1, 2 ,...: 


Ps{tN) 



(60) 


Here we have introduced the tensor S, 


W. 



(61) 


where Un'^ is shorthand for Un+(i-i)M and V® for V„_|_(;_ 2 )m- 

We can think of as a superoperator acting on k copies of the system which we can label Si , 
for Z = 1,..., fc. A compact way to write psi^N) can be found by introducing a time-propagator 
on a Zc-fold system spaee, ^m(Z 7 v): 


Si - 

SM{tN) = : 

Sk-i — 

' ' 

p(k—l) 

- 

— ■■■ — 

• 

p(fc —1) 
^N-kM + 1 


' 

ci^) 

^N-kM 

— ■■■ ~ 

^(fc) 

Ok 



n{k—l) n{k—l) ci^) ci^) 

<'M • • • ^N-kM+l^N-kM • • • 


( 62 ) 
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Ps(^Af) is found by imposing the following boundary conditions: 

Ps(ijv) = .. .Tr(S2.5i)(63) 

where the generalized trace, Tr, is defined in Eq. (22), and ps^ is the system initial state. 


B. The continuum limit 


We wish to take the continuum limit of Eq. (62). Again, we assume for simplicity an incoming 

(k'\ 

vacuum field, with cr„ = |0)(0| for all n. We first note that Sn can be built from the following 


tensors: 


? 




/" N 


u, 


(fc) 




= tr„(fc)(•®o-n(fc)) 

r{A;)i 


= exp I ( 


(64) 


(65) 


where we have defined n(l) = n + (Z — 1)M and the superscript denotes the system on which an 
operator acts. 


( •) =- - • = trn(O) ® cr„(0)) 


= exp • +MV[l\^'^]> 


( 66 ) 

(67) 


and the composite tensor 




7 ;d+l,0(.) = 


- 


I 




( 68 ) 


V, 


(*+ 1 ) 




= exp (^ - i^‘H[H''g+^'> + Hf] • +MV[Lf>] • +AtT'[4'+^^] • +AzC[l('\ 4'+^’] • j, (69) 

where 

c[x,y] • = A . yi + y • At - ytA • - • ytA. 


(70) 


can be written in terms of these tensors as a product 

f(k) _q-{k)q-{k,k-l) q-{3,2)q-{2,l)q-{l) 

'-'n. ' n. ' n ' ' ’ ' n ' n ' n 




where, for notational convenience, we have defined 

+ 77^+^’ + i(Lt'tti,6+i) _ H.C.), 
L,,,+i = L«+4'+'>, 


(71) 


(72) 

(73) 
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except for 


Ho,i 




{k) 
S ’ 


T _r (1) 

-^ 0,1 —-^2 ’ 

T _r (^) 

J^k,k+l — -^1 

We can now write the continuum limit of Eq. (62) as 


E^{t) = lim 
N—yoo 


(k—i) p(fc—1) p(fc) cy 

• • • ^N-kM+l^N-kM •••<-! 


(k) 


1=0 


=rexp( [ ds^\-lH[Ht^\s)+H^i\s)]+V[U,i+,{s)] 


where T is the time-ordering operator, and we have defined 

^W(s)=ylW foralU^fc, 
ylW(s) =0[t- (fc- 1 )t-s]^W, 

for any system operator A, where 6{t) is the Heaviside step function. 
We can differentiate Eq. (78) to find a differential equation for Sa{t): 


ds 


1=0 


= E -9^ + V [L,,(+i(s)] f,(f). 


(74) 

(75) 

(76) 

(77) 


(78) 


(79) 

(80) 


(81) 


with initial condition fo(t) = 

The continuum solution for ps{t) is found from £r{t) by imposing the boundary conditions in 
Eq. (63): 


Ps{t) - ^J:{Sk,Sk-i) •••Tr(S2_Sj)5^(t)pSi, 


(82) 


where ps^ is the initial state. Eq. (81) and Eq. (82) are the main results stated in the letter. 
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